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ABSTRACT 

In this article we revise the problem of anomalous values of pulsars' braking in- 
dices Uobs and frequency second derivatives iy arising in observations. The intrinsic 
evolutionary braking is buried deep under superimposed irregular processes, that pre- 
vent direct estimations of its parameters for the majority of pulsars. We re-analyze the 
distribution of "ordinary" radio pulsars ona iy — v, iy — v, ii — v and nob^ —TcK diagrams 
assuming their spin-down to be the superposition of a "true" monotonous term and a 
symmetric oscillatory term. We demonstrate that their effects may be clearly separated 
using simple ad hoc arguments. Using maximum likelihood estimator we derive the 
parameters of both components. We find characteristic timescales of such oscillations 
to be of the order of 10^ — 10'* years, while its amplitudes are large enough to modulate 
the observed spin-down rate up to 0.5-5 times and completely dominate the second 
frequency derivatives. On the other hand, pulsars' secular evolution is consistent with 
classical magnetodipolar model with braking index n 3 

So, observed pulsars' characteristic ages (and similar estimators that depend on 
the observed v) are also affected by long term cyclic process and differ up to 0.5-5 
times from their monotonous values. This fact naturally resolves the discrepancy of 
characteristic and independently estimated physical ages of several objects, as well as 
explains very large, up to 10® years, characteristic ages of some pulsars. 

We discuss the possible physical connection of long term oscillation with a complex 
neutron star rotation relative its magnetic axis due to influence of the near-field part 
of magnetodipolar torque. 

Key words: stars: pulsars: general - methods: statistical 
PACS 04.40.Dg - 95.75.Pq - 97.60.Gb - 98.62.Ve 



1 INTRODUCTION 

Radio pulsars are highly periodic variable astrophysical 
objects, powered by neutron stars' rotation with periods 
evolving with time. Observational determination of their 
timing parameters - instantaneous p eriod and its deriva- 
tives - is an extremely complex task (|Edwards et al.l [20061 : 
ICordes fc Shan non 2010J, but after the correction of pulse 
times of arrivals (TO As) for the effects of radio wave propa- 
gation in the interstellar medium, motion and position of the 
Earth, gravitational delays in Solar system potential etc., 
the resulting phase of the light curve may be well described 
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by an (infinite) Taylor series 

<i>{t) = <^o + v(t - to) + ii>(t - to)^ 4- i/>(t - to)^ + ... (1) 

dominated by the lower order terms. This may also be ex- 
pressed in terms of observed rotational frequency as 

v(t) + v{t - to) + h>{t - tof -f i ■/>■(* - to)^ + ... (2) 

where values of v and its derivatives are attributed to the 
epoch to. 

At the same time, the majority of published theoretical 
models of pulsar braking predict the spin-down described 
by a differential equation of the form (jManchester fc TavloJ 
ll977l : lBeskin et al.lll993l i 

V = -Kv"^, (3) 

where n is a (constant) braking index and K depends on the 
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individual physical properties of the pulsar. For the vacuum 
magnetodipolar model the canonical value is n = 3, pulsar 
wind decreases it to n = 1 and mu ltipolc magnetic field 
increases it to n J5 5 ^Manchester fc Ta ylor f 977i ). 

For such a spin-down with constant K the braking index 
is equal to a simple combination of frequency and two its 
derivatives Q 



vv 

nobs = ^ = n 



(4) 



That is why the determination of spin frequency second 
derivative via pulsar timing is important for understanding 
the pulsar' spin-down. 

Several decades of very de tailed timing of hu ndreds 
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julsars (D'Alcssandro et al. 1993; Bavka l et al.l 119991: 



Chukwudc 2003 : Hobbs et al...2004 . .201Q, : Tivingstone et all 
2005^■ including the best studied one - Crab pulsar 
I Scott et al.ll2003 ) - demonstrated, however, that the mea- 



sured rotational phase ((!]) evolution is not generally consis- 
tent with the one expected for a braking law @ with phys- 
ically reasonable parameters. Observed rotational phase (j) 
(as well as v, v) includes components with a very complex, 
irregular behaviour. The first kind of such irregularities is 
"glitches" - sporadic fast changes of pulsars' periods and 
spin- down rate with up to f ew mo nths relaxation timescale 
(e.g. [Manchester fc Tavloil (ll977D ). while the other one - 
quasi-random phase variations with typically red Fourier 
spectrum and characteristic timescales of a few months to 
few years - "timing noise ' ' (e.g Cordes fc Hclfand (1980); 
iLvnd (|l999l ): lHobbs et all (|2010l )). 

Both of them affect the measurements of i^, and u, and 
often make them dependent on the duration and epoch to of 
the observations' time span, when it is shorter or compara- 
ble to timescales of these processes. But, if observations are 
performed over sufficiently long intervals of about a decade 
or longer, the observed coefficients of series ([1} become more 
stable - the influence of these irregularities is mostly sup- 
pressed (see Sec. 12.11 for additional discussion). To date, 
there are more than 200 pulsars observed fo r longer than 
15 years llfiavkal et al.lll999l : IChukwudd [20031 : iHobbs et al] 
12004 |2010| ). and still the majority of i> values, measured 
over such long intervals, lead to extremely high, up to 10^ 
(and even more), braking indices. Moreover, their values are 
negative for nearly a half of all objects. In general, measured 
braking indices drastically differ from physically reasonable 
values, and this raises serious problem on the pulsars' long 
term spin evolution. 

While glitches are widely believed to be caused by a 
rapi d transfer of momentu m from the neutron stars interior 
fe.g. lEspinoza et al.l (|201ll ) and references therein), the na- 
ture of the timing noise is still unclear and no self-consistent 
and widely accepted model of pulsars' phase residuals irreg- 
ularities is constructed. There have been numerous attempts 
to attribute it to various s tochastic phenome n a in the neu- 
tron star magnetosphe re (jCheng et al.l Il987 Contopoulos 
l2007l : iLvne et al.1 12010| ). in the interior (|Cordes fc Dowiis 



^ Note that the observed ones, estimated by fitting the TOAs 
with a model of JTJ broken down to the third order term, are 
always biased even if the spin-down is exactly according to equa- 
tion ^ due to influence of non-zero higher-order terms. For actual 
pulsars, however, this bias is negligible and will be ignored below. 



Il985h . to spin-down torque variations (jUrama et aLlliooil ). 
and to the existence of a number of differe nt spin-down 
regimes for a single pulsar (|Lvne et al.ll201(il ). Purely phe- 
nomenological attempts to describe observed noise as ran- 
dom walks of different orders (in pha se, frequency or its 
derivative) (|Cordes fc GreensteinI [l98ll i have also not suc- 
ceeded as with the increase of observations time spans 
the noise appeared to be more complex than these simple 
models predict. Some characteristics of this essentially ir- 
regular process, however, have been parametrized through 
sev eral noise strength parameters and their correlations. 



So, [a rzoumanian et al] (|l994f ) parametrized pulsars' timing 
noise through Ag oc logdi^l/z/) quantity and found a good 

Ag — P correl ation, where P is a pulsar p eriod. 

In turn, ICordes fc Helfandl (|l980l ). ICordes fc DownsI 

l|l985h and then Chukwudd ( 200; ) have investigated activity 
parameter A, which is a direct measure of timing residuals 
variance. A good correlations of A versus !> and D was also 
found. 

However, the quantities used in these papers to 
parametrize timing noise were in fact just the measures of 
observed frequency second derivative. Even the activity pa- 
rameter A was calculated using timing residuals after sub- 
traction of a second order polynomial in series ([l|. There- 
fore, the correlations obtained simply reflect the quite strong 
D — u dependence. 

At the same time, only a few long-term mechanisms 
have be en offered to explain obs e rved anomalous braking 
indices. iDemiariskv fc Proszvhskil (|l979l ) studied the varia- 
tions of pulsar timing parameters due to existence of pos- 
sible massive partner - most of radio pulsars are, ho wever, 
believed to be isolated objects. I Alpar fc Bavkall l|2006l ) inves- 
tigated the impact of unresolved or missed glitches on the 
observed timing parameters, but it is very difficult to explain 
some extre mely hig h braking indices \n[> 10^ — 10^ by such 
mechanism. iGullahorn fc RankinI (ll977|) considered varia- 
tions of pulsar braking torque with timescales of 10^ — 10'' 
years, probably due to interaction of the pulsar with inter- 
stellar medium in its vicinity. The idea of spin-down torque 
variations itself is quite promising and reasonable, but the 
nature of variations on such timescale was not pointed out 
in their note and remained un cl ear. 

In turn, in iBeskin et al. 1 l|2006l ) and iBirvukov et al.l 
( 2007) we considered the existence of a long-term cyclic vari- 
ational process affecting pulsars' spin-down on a timescale of 
thousands of years. One may ca ll such a process a thi r d type 
of timing irregularities. Later, iBarsukov fc Tsveanl (|2010D 
suggested that the physical nature of such a process may be 
related to the "anomalous" magnetodipolar braking torque, 
which may produce a forced precession of pulsar's rotational 
axis around its magnetic moment with the period of 10"^ — 10^ 
years, and demonstrated the possibility of significant varia- 
tions of obse rved rota t ional parameters on such timescale. 
Even earlier, iMelatosI (|2000l ) also investigated the triaxial 
neutron star rotation taking into account this torque and 
demonstrated that the rotation of such star may be very 
complex and even induce variations of magnetic inclination 
angle, which will also affect the neutron star spin-down. 

As a further development of this concept of a long-term 
variational process, in the current work we analyse the en- 
semble of 297 pulsars with published data on second fre- 
quency derivatives and demonstrate that these values may 
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be used to estimate the parameters of both secular (evo- 
lutionary, monotonous) and additional, cyclic, components 
of a pulsar spin-down under simple and reasonable assump- 
tions. We formulate such a two-component model of pulsar 
braking and derive its parameters using a maximum like- 
lihood estimator. The results are quite reasonable, as the 
parameters of monotonous component are in a good agree- 
ment with a standard magnetodipolar spin-down model. 

The article is organized as follows. In Section [2] we jus- 
tify that measured i> values can indeed be used to analyse 
pulsars' spin-down. In Section [S] statistical analysis of a 297 
objects is performed, phenomenology of a two-component 
pulsars' spin-down is presented and parameters of its irreg- 
ular term are roughly estimated using model-independent 
arguments. In Section |3] this two-component model is used 
for the determination of the secular spin-down parameters 
also. In Section[5]we discuss the results and its astrophysical 
implications. The conclusions are given in Section [S] 



2 PULSARS' i> MEASUREMENTS 

2.1 What i/'s can tell us on the pulsars' physics? 

Anomalous values of /> (and braking indices) for most of or- 
dinary pulsars raise at least two principal questions, critical 
for any serious analysis of these quantities. Namely: 

• Are values measured over long timespans (decades) re- 
ally related to pulsars' or interstellar mediuir|j physics or 
they are just artifacts of somewhat incorrect observations 
and/or data reduction? 

and 

• If they are not artifacts, then is the anomaly induced 
by the same phenomena responsible for the short time scale 
(months) quasi-periodic timing irregularities with red-like 
spectra, observed in timing residuals? 

The answer to the first question above is clear. Pul- 
sar timing is a complex, but well-defined and thoroughly 
described procedure, and the sources of observational un- 
certainties arising in it are well known and are taken into 
account during the reduction, as well as various effects that 
affect the timing systematically - e.g. motio n of the Earth, 
relativ istic effects, etc (see, for example, lEdwards et al.l 
l|2006t )'). There is no reason to believe the significant cu- 
bic tren ds, clearly seen in phase solutions for hundreds of 
pulsars (|Hobbs et al.ll20o3 ). are artifacts. 

There is also a method for indirect measurements of 
a pulsar' second derivative by comparing !> values acquired 
over sufficiently distant, sepa rated by 6—20 years, time spans 
ij Johnston fc Gallowavl[l999l ) . Results of such estimation are 
anomalous too, and are generally consistent with directly 
measured values. 

Obviously, all accurately measured i> values, even being 
anomalous, are physically meaningful - they are indeed due 



The effects of propagation of pulsar radiation through its "near" 
and "far" medium can introduce an additional delays to the 
times of pulses arriva l s and therefore affect the timing solution 
dCordes fc Lazidl2002l : [ Cordes fc ShannonlbOlOl) . 



to peculiarities of pulsars' spin down, while their anoma- 
lousness is just an indication of insufficiency of our models 
of pulsars' braking. 

On the other hand, the answer to the second question 
is not so obvious, as stochastic timing irregularities, seen in 
a large fraction of pulsars, still are not well understood and 
are complex phenomena. 

Phenomenologically, timing residuals' irregularities and 
unexpected anomalous parameters of timing solution have to 
be separated in the analysis. Indeed, the former ones are seen 
directly within observational datasets and their stochastic 
nature is clear while the origin of the latter ones can only 
be speculated about. Below we argue in favour of indeed 
different origins of these phenomena. 

The phase residuals after all appropriate corrections 
and quadratic trend subtraction (i.e. after extraction of f 
and !> only) form a very plural zoo over the pulsar popu- 
lation. There are at least several diff erent classes of typ ical 
pulsars' timing series (see Figure 3 in lHobbs et all l|2010l )') - 
the ones with dominant influence of a cubic term (e.g. PSRs 
B0959-54, B0943-I-10, B1657-13 etc.), the ones with more 
complex but generally smooth behaviour (PSRs B1620-26, 
B1706-16, B1727-47 etc.), and the ones that are purely white 
noise supposedly due to measurement errors. The latter class 
typically contains pulsars with inaccurate or unmeasured 
second frequency derivatives, and will not be analysed here. 
The second one is often used to depict the timing deviations 
from an expected spin-down law as a process with red noise- 
like spectrum; there are, however, some evidences that this 
behaviour in fact is a combin ation of a cubic ord er term and 
a quasi-periodic component (|Hobbs et al.|[201ol ). Moreover, 
young pulsars (e.g. Crab, B1509-58, B2011-I-38) show strong 
short-term quasi-periodic timing noise, along with cubic 
trends corresponding to a spi n-down with phys i cally mean- 
ingfu l braking indices ^ 3 I Lvne et al. _ 19931: IScott et al.l 



I2OO3I : [Livingstone et al.ll2005l : iHobbs et al.ll20ld ) 



Also, the behaviour of i>'s with the increase of observa- 
tional span is not consistent with their origin due to stochas- 
tic or short timescale noise processes. For example, for PSR 
B1706-16, variations of i> with an amplitude of 10~^* s~^ 
have bee n detected on a tim escale of several years (see Fig- 
ure 7 in iHobbs et all (|2004l )') - its values, however, are al- 
ways around the one revealed by the fit over the entire 25 



year time span (i> — 3.8 x 10 
braking index ^ 2.7 x 10^. 



s ), which still leads to a 



All this strongly suggest that timing noise is a process 
distinct from the one producing large (and anomalous) cubic 
trends in ((T} ; this noise acts mostly as a randomizing factor 
decreasing the accuracy of i> measurements, but not defining 
their properties in general. 

As a result, we believe the second derivative values to 
be very promising tool for studying long timescale features 
of pulsars' spin-down. Although the i> values of individual 
pulsars bring little information on the spin-down of neutron 
stars, the properties of its distribution over an ensemble of all 
pulsars seems to be appropriate for spin-down study under 
a simple and reasonable assumptions on the properties of 
process producing these anomalous cubic trends. 
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2.2 The subset 

The set of pulsars under investiga tion is simil ar to 

the one used in ou r previous works (| Beskin et alj l200d : 
iBirvukov erahl 120071). From the 39 3 objects of the ATNF0 
catalogue ( Manchester et al]|2005l ) with known u we have 
compiled a list of "ordinary" radio pulsars that 

• have P > 20 ms and jPj > 10~^^ s/s (i.e. had not been 
recycled) ; 

• have relative accuracy of second derivative measure- 
ments better than 75%; 

• are not components of binary systems and not in a list 
of known anomalous x-ray pulsars. 

19 supplementary pu ls ars from oth er sources 
l|D'Alessandro et alj 1 19931 : IChukwudd |2003| ) have been 
added . The final set con sists of 297 objects including 247 
from l|Hobbs et al.l [20o3 ') . 18 of them are associated with 
young supernova remnant^ 



3 STATISTICAL ANALYSIS OF PULSARS' 
TIMING PARAMETERS 

3.1 V — V — V dependency. Cyclic evolution of 
pulsars' spin-down? 

Our work is based on the study of relations between quan- 
tities derived from pulsars' timing: spin frequency, its two 
derivatives v and observed braking index Uohs = vv/i/^ 

and characteristic age t^h = —v/2i'. 

As we previously demonstrated in l|Birvukov et al.l 

I2OO7I ). the logarithms of li^l show significant correlation with 
the ones of — /> for both positive (172 pulsars, correlation 
coefficient r ~ 0.9) and negative (125 pulsars, r ~ 0.82) 
branches of the v — v diagram (Figure [1} . 

The apparent separation of branches on the v — v dia- 
gram is primarily due to the logarithmic scale of the plot, 
and actually objects continuously cover the full range of v 
values, excluding only the gap on small ones due to limited 
accuracy of measurements (which is not better than 10~^^ 

Young pulsars confidently associated with supernova 
remnants are systematically shifted to the left on the di- 
agram (open symbols in Figure [T]) and are absent on the 
right. Hence, pulsars seem to evolve towards lower values 
of \v\. Moreover, if t^k is an appropriate pulsar age estima- 
tor, then \v\ is the same too, as they are highly correlated 
(r = 0.99 in logarithmic scale, see Fig. O 

On the V — V diagram, shown in Figure O pulsars' be- 
haviour is similar: they are born with higher values of v and, 
since z> < 0, evolve toward lower values. The direction of the 
evolution is the same for positive and negative branches of 

V — V. 'SiQ all older pulsars have systematically lower values 
of which is consistent with evolutionary interpretation of 

V — V diagram. 

So, we conclude that each pulsar during its evolution 
moves along the branches of the v — ii and v — v diagrams 



^ http: / /www. atnf . csir o ■ au/r6search/pulsar/psrcat/[ revi- 
sion from 6th Oct 2007 

* Data on PSR-SNR associations was also taken from ATNF. 
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Figure 5. The correlation between i/ and characteristic age 
= —vjiv for 297 pulsars with measured />, as well as for 
1337 ordinary isolated pulsars taken from ATNF database. The 
correlation coefficient between logarithms of these quantities is 
r = 0.99. The spread of points on the scatter plot is due to the 
distribution of frequencies u among the pulsars. It, however, is 
unable to break tight relation between t^^ and u. Therefore, the 
frequency first derivative may also be considered a characteristic 
of pulsar age. Additionally, the obtained dependency —u oc r°,^ 
with a 7^ — 1 suggests the presence of significant i> — f correlation 
(see Fig. |4|| 



while its !> value increases and v value decreases. However, 
on the negative branch of iy — v, the value of the first deriva- 
tive, being negative too, can only decrease with time, since 
is a formal derivative of z>, and both of them are regular 
components of the observed rotational phase ([T}. So, pul- 
sar motion along the branch may only be backward, which 
clearly contradicts its evolutionary interpretation suggested 
earlier. The solution we offer is to assume a cyclic behaviour 
of pulsars on this diagram. As pulsars evolve, they repeat- 
edly change sign of i>, in a spiral-like motion from branch to 
branch, and spend roughly half their lifetime on each one. 
Such behaviour is sketched in Figure (2] 

The timescale of such variations has to be much shorter 
than the pulsar life time, and at the same time significantly 
larger than the one of the observations to systematically 
affect the timing solution. Of course, such cyclic behaviour 
should also affect other spin-down parameters - f and v. 

Then, on v — v plot, shown in Figure |4l pulsars move 
from bottom right to the upper left, from high to low f and 
|!>|, forming quasi-evolutionary sequence, with linear regres- 
sion coefficient ~ 2 in logarithmic scale. Unfortunately, this 
diagram itself can't be used to study the spin-down evolu- 
tion of pulsars, as it is distorted by various selection effects 
(death line crossi ng, emission beam w idth dependency on 
pulsar frequency (jTauris fc Manchester! 1998'). etc) and cor- 
relation of pulsars' initial v and v at birth. 



3.2 Phenomenology of observed pulsars' 
spin-down 

Now we introduce an appropriate phenomenological descrip- 
tion of the complex timing evolution of pulsars including the 
cyclic variations we proposed above. 

The exact form of a pulsar's trajectory on the iy — v dia- 
gram is unknown, but in general its rotational evolution may 
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Figure 1. The u — u diagram for 297 pulsars. The figure shows the pulsars taken from lHobbs et al" as circles, and the objects 

measured by other groups as squares. Open symbols represent relatively young pulsars confidently associated with supernova remnants. 
Analytical fits for both positive and negative branches are shown as solid lines. They were obtained by means of linear least-squares 
regression in the logarithmic scale. Measurement errors are shown as error bars and in most cases are well inside the symbols. The 
diagram is an evolutionary sequence - pulsars systematically move from left to right of it. Two branches are due to cyc lic variations of 
pulsa rs' rotational parameters on a timescale significantly shorter than lifetime but longer than time spans of observations (iBirvu kov et al.l 
[2003) as illustrated in Figure [2] The dashed lines represent the power-law approximations of the upper envelopes of positive and negative 
branches, estimated by the method described in IC^dielliooi). Their moduli may be used as an upper and lower limits on the variational 
amplitude of u, respectively (see Section l3.4p . 
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Figure 2. Concept illustration of a pulsar motion on the u — v diagram due to a cyclic evolution of the timing parameters. Monotonous 



component of a spin-down (dashed line) is according to a canonical law i>ev = ~Kv^^ with K = 
Superimposed cyclic spin-down component is of a 5!>(t) = v^^ cos(f!i) form, with A = 0.2 and 27r/f7 



10- 



and i'eti(O) 



30 Hz. 



10^ years. Resultant track of a 



pulsar is shown as the solid line. On the left panel, an initial stages of pulsar evolution are presented - the track is shown for 10 — 10 s 
interval of ages. Here the amplitude of v variations is not much different from evolutionary value v^v, and therefore the youngest pulsars 
are unable to change the sign of their observed u. On the right panel the evolution of the same pulsar during (1 — 1.3) X 10^^ s ages 
interval is shown. Here, v variations prevail, and the pulsar repeatedly changes the sign of i). Note, that due to positive contribution of 
variations are never exactly symmetric in respect to v = 0. 
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Figure 3. i) — v diagram for the subset of pulsars under investigation. This diagram is generally similar to the v — ii diagram shown in 
Fig. [T] and has the same evolutionary meaning, as pulsars born with higher values of v and evolve to the lower ones. Dotted line shows 
typical evolutionary trend for a pulsar evolving according to !> = —Ku'^ law with K = 2 X 10~^^ and n = 3. 
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-2.82 X 10" 



is also shown according to l|Bhattacharva et ahlFigQah 
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be described as the superposition of the regular component 
of the observed rotational phase (j>evit), and an irregular, 
cyclic one S(j){t): 

4>{t) = 4>,^{t) + 5(j>{t), (5) 

where 

/ 5(t>{t)dt « (6) 

J AT 

over a long enough time interval AT comparable to a pul- 
sar's lifetime. After term by term differentiation of equation 
([SJ one gets: 

u{t) = u^^{t) + 5jy{t) (7) 
for the rotational frequency and 

U{t) = + Sv{t) = Uev{t) [1 + e(t)] (8) 

for the spin-down rate v. Obviously, due to secular loss of 
rotational energy, the value of v^v {t) should be always neg- 
ative for an isolated pulsar. Quantity e{t) = 5i'{t) /uevit) is 
a spin-down rate relative deviation, which is not necessarily 
small, and is likely greater than -1: 

eit) > -1 (9) 

due to the absence of "ordinary" pulsars with !> > in the 
subset under investigation. It is clear that deviation e(t) is 
also cyclic; in the simplest case it is a harmonic function of 
variational phase ifiit): 

e{t) = A cos <p{t), (10) 

where ifi{t) = (po + Qt and A is the relative amplitude of the 
!> oscillations, related to the absolute amplitude Av'- 

Ai,{t) = Ai'evit) (11) 

By second differentiation of equation ([5]) for a rotational 
frequency second derivative we get: 

U(t) = Uev{t) + Suit) = l'ev{t)[l + Tj{t)] (12) 

where 2> Vevit) (i.e. \'r]{t)\ 2> 1) for most pulsars, 

and this is the reason behind anomalously high values of 
observed \D\ and braking indices. In other words, the abso- 
lute amplitude of the i) variations An is much greater than 
the regular term Dev Note also that cyclic terms e(t) and 
rjit) axe not independent - they are connected through the 
relation: 

r^{t) = e(t) + mp^ (13) 

Since pulsars secularly evolve to the higher values of i> 
and i>ev{t) < 0, the sign of Devit) ~ di'ev{t)/dt should always 
be positive: 

l'e.v{t)>0 (14) 

Such positive contribution should introduce a small 
asymmetry of the observed values of D in respect to u = G, 
even if the variations of Sv axe intrinsically symmetric. Such 
asymmetry drives the average motion to the right on the 
D — ii diagram and affects the times that pulsars spend with 
positive and negative D. 



3.3 Asymmetry in observed i^'s 

Indeed, numbers of pulsars with positive (A''^ — 172) and 
negative {N~ = 125) values of D within the subset are sig- 
nificantly different. If this difference is just accidental, its 
probability is too small, P = 6.4 x 10~^ only, assuming bi- 
nomial distribution of 297 objects with p = 1/2 chance to 
be on a positive branch. Therefore, null hypothesis of ex- 
actly symmetric branches is rejected on a 0.64% significance 
level, and the branches are indeed asymmetric in number of 
object 

We investigated in detail the behaviour of such signifi- 
cance levels for a number of subsets of pulsars with v greater 
than and v less then a given value. Corresponding dependen- 
cies are plotted in the Fig. (6] It is clearly seen that branches 
of the V — V and v — v diagrams are significantly different 
only when relatively young pulsars are taken into account, 
while behaviour of the older pulsars seems to be the same 
for both signs of v. It is consistent with the existence of a 
positive evolutionary trend v^^, which is negligible for older 
pulsars, but is large enough to affect the second derivatives 
of the younger ones. 

In other words, for the initial stages of pulsars life, val- 
ues of \hv\ are less than i^eu (i.e. |??(f)| < 1), and the second 
derivatives are always positive (see the upper left region of 
Fig. [!}. Later on, as \v\ decreases, the v^^ becomes quite 
small and the relative deviations \r\(t')\ grow up extremely 
over unity; as a result, observed values of v turn out to 
be significantly different from an intrinsic v^^, and it corre- 
sponds to the great increase of their absolute braking indices 
(see Fig. O. These different stages of pulsars' evolution are 
illustrated in the Fig. [2] 

At the same time, the absence of a significant difference 
between the v — v branches for relatively old pulsars - their 
symmetry - implies an important fact that variations are 
indeed approximately symmetric in respect to evolutionary 
trend: a pulsar deviates to higher and lower values of v with 
nearly the same amplitude and spends an equal amounts of 
time with v greater and less than v^^. This fact will be used 
in Section m for the estimation of the parameters of pulsars' 
two-component spin-down model. 

However, some generic characteristics of the cyclic, ir- 
regular component of v behaviour - its amplitude and 
timescale - may be easily derived from an exploratory anal- 
ysis of the V — ii and Uoba — TcK scatter plots using simple, 
model-independent arguments. Such analysis is presented in 
the next two subsections. 

3.4 Simple limits on the timescales and 
amplitudes of the cyclic process 

As stated above, for the majority of pulsars the amplitude 
Ai, of the second derivative variations is much greater than 
the positive evolutionary trend Ve.v AH the pulsars on the 
negative branch of i> — ;> diagram are due to subtraction 
of this amplitude from the trend. Therefore, the values of 
\v\ of the lower envelope of this branch (in terms of \v\) are 
sensible estimation of a lower limit on At, . At the same time, 

^ At the same time, non-parametric Kolmogorov-Smirnov test 
rej ects the hypothesis of the common distribution of | i> | of positive 
and negative branches on a 2.5% significance level. 
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Figure 6. The significance levels of a non-parametric binomial 
test for the hypothesis of equal numbers of pulsars with u > 
and u < 0. Each point on the plots represents the p- values for the 
subset of pulsars with i> greater than (upper) and i/ less than (bot- 
tom) selected threshold. Minimal significance levels correspond to 
the complete set of pulsar. So, the differences of branches in num- 
ber of pulsars are due to contribution of the youngest objects. The 
absence of any significant branch differences for older pulsars fa- 
vors the assumption of symmetric variations of i> in respect to its 
evolutionary values. 



the values of D of pulsars near upper envelope of positive 
branch may be considered an upper limit of An . 

Indeed, due to existence of the secular trend I'ev > 0, 
the amplitude of the variational term Si> can not exceed the 
values of i>'s of the upper envelope of the positive branch (for 
a given !>). For pulsars there i> = +Si} = D^v + An > Av, 
while for ones on the lower envelope of the negative branch 

\i>\ < Ai> — Vev < Ai>. 

We used the Bou n dFit generalized least-squares code 
developed by ICardiell (|2009l ) to acquire an analytical ex- 
pressions for the mentioned envelopes of the negative and 
positive branches on the v — v diagram, and found them to 
be 



and 



10 



10^ 



-9.51/ . nI.05 



(15) 



(16) 



respectively. They are shown as dashed lines in Figure [T] 
Since we consider pulsars on fits (|15p and p6|) to be ob- 



served with amplitude values of frequency second derivative, 
their i/'s are close to turning points, where corresponding (5z> 
change signs during the oscillations (see also the sketch in 
Figure [2]), and are obviously small (iv ~ 0). Hence, the as- 
sumption of ~ !>5i, in equations (|15[) and (|16p is justified. 

Due to the absence of pulsars with positive the 
amplitudes of variations of frequency first derivatives are 
likely to be less than their secular values: Ai, < —Um and 
Ai,, up ~ —Vev Also, for any more or less stable cyclic pro- 
cess, the following relations should be true: 



(17) 



where = 2-n /T is a characteristic frequency of cyclic vari- 
ations. (For the special case pop these relations are exact). 
For the upper limit Q,up ~ An ,up / Ac ^lom while for the lower 
one Q,iow ~ Ai)^iow/ Av^up- Hence, one can estimate: 

27r 



i Hon 



< -27r- 



Ai} ^low 

which leads to upper limit on the timescale Tup'- 



Tu 



(1.2 X 10*^ years) • (-^>e^,,l4)°■°^ 



(18) 



(19) 



where i/^^ is normalized to 10 s ^. 



The value of T,, 



changes from ~ 8 x 10^ years for oldest to ~ 2 x 10® years for 
youngest pulsars. Then, the value of lower limit Tiow should 
depend on the upper limit of An and typical minimal relative 
Sv amplitude Am. Ac,, low = —AmVev Hence 



Tlou 



-27r 



Tlori 



(650 years) • ■ (— !>ei.,i4)" 



(20) 



(21) 



However, T should be more than few times of typical ob- 
servational span length for the pulsars under investigation, 
which is 15 — 20 years. Hence, 



Tlou, ~ 50 - 100 years 



(22) 



and Am is unlikely less than ~ 0.1. Moreover, we will show 
below that Am is even greater than ~ 0.5. 

At the same time, assuming that fit to negative branch 
of i> — ;> diagram 



-10" 



-10" 



(23) 



represents in the same way a typical variational amplitude 
(shown as solid line i^-(i') in Figure[l|, one may, in a similar 
way, estimate a typical timescale able to provide observed i> 
spread. So, Qtyp ~ \v-\li'ev or 



Tty 



2tt 



n 



typ 



(7.5 X 10'' years) ■ (-i>ei,,i4)°-°^ (24) 
9 X lO** years for oldest and 



with spread from ~ 6 x 10^ to ' 
youngest objects, respectively. 

Finally, according to equation (I17p the amplitudes of the 
pulsars' frequency variations should be roughly times 
the ones of its derivative and times - its second one. 
For typical maximal timescale ((24]): 

A, < v1j\V- \ = (4 X 10-=* Hz) ■ (-!>e„,14)"'''^ (25) 

These values do not exceed a few Hertz even for young pul- 
sars and significantly less than observed frequencies for most 
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of pulsars. It justifies the neglection of variational term in 
the rotational frequency decomposition ((TJ: 

(5l/(t) < T^evit) and u{t) ~ Vevit) (26) 

3.5 Pulsar observed braking indices and 
characteristic ages. An analysis of the 
nobs - Teh diagram. 

There are number of widely accepted estimators of pulsars' 
physical properties based on results of timing solution and 
power-law ((S)) spin-down model with constant K and n — 3: 
characteristic age Tch = —i'/2i', rotational energy loss rate 
Erot oc vv and surface magnetic field 

B = (3.2 X lO"-* G)\/-!>i^-3. (27) 

But, if the amplitude of frequency derivative variations, 
|e(f)|, is large, these quantities may be significantly biased, 
as they are computed using observed v values, and not actual 
evolutionary ones Oev 

Using phenomenology introduced earlier and assuming 
power-law spin-down defined by equation (|3]) with constant 
K and evolutionary braking index n, observed quantities 
Uobs and Tch may be expressed as: 



and 

r..= (^^t+^j^ = ^, (29) 

where Po is pulsar initial period, t - its real age, and e and 
rj are relative i> and i> deviations according to equations ((8]) 
and mi). If n = 3, then /2K ~ 10^ years for the typical 
pulsar and r^h.ev ~ i for the older objects. 

So, if e varies in ±0.9 range, then half of pulsars will 
be observed as up to Tch/Tch,ev ~ f/(l ~ 0.9) = 10 times 
older than their true unbiased characteristic ages, while an- 
other half as 1/(1 -I- 0.9) « 0.5 times younger. The effect is 
therefore seemingly asymmetric, as characteristic ages will 
mostly appear overestimated. 

The same bias will also appear in the rotational energy 
loss rate, while magnetic fields will be 2 — 3 times underes- 
timated for pulsars with e < 0. 

Moreover, since observed braking indices Uoba = 
vuv~^ oc z>~^, their absolute values will be additionally in- 
creased up to order of magnitude for a half of pulsars. 

Since Uobs and Tch depend on e(t) and rj{t), the pul- 
sar behaviour on the Uoba — Tch plane, shown in Figure [3 
is similar to that in the D ~ u one. The positive and nega- 
tive branches are not identical: the numbers of pulsars and 
shapes of branches are significantly different. 

The most interesting regions of the diagram are the ar- 
eas I, II and III with the relatively young pulsars. There are 
19 pulsars with positive and only 2 with negative riobs in the 
area Ilfl Binomial test with p = 1/2 rejects the hypothesis 



of an accidental origin of such asymmetry on a 0.01% signif- 
icance level. The pulsars with positive Uobs have measured 
Wobs ~ 10 — 50, which are larger than any reasonable sec- 
ular value. Moreover, they slightly deviate from the rest of 
positive branch of the diagram. 

Is it possible to explain these a result of decay 

of K coefficient in spin-down law ((3}, pr obably due to evo- 
lution of magnetic inclination angle (e.g. lDavis fc GoldsteinI 
(|l970l ))? To yield pulsars with both Uobs ~ 30 - 50 and 
Tch ~ 10* years, corresponding timescale —K/K should 
be as short as few thousands years only. However, e.g. 
iFaucher-Giguere fc KaspH (|2006| ) have shown that there is 
no significant decay of coefficient K on a timescale ~ 10* 
years for isolated pulsars. 

We plotted pulsar evolutionary tracks corresponding to 
spin-down law ^ with n = 3 and —K/K = 5 x 10* and 
10^ years on the Uobs — Tch diagram (blue dashed lines in 
Figure (Tjl. It is clearly seen that first of them is insufficient 
to explain Uoba of pulsars in area III, while second one is 
unable to affect observed distribution of pulsars' Uobs at all. 

So, we conclude that such asymmetry results from com- 
bination of a small enough amplitudes of u variations (i.e. 

<C 1, so the negative values of i> can not be reached yet) 
and high amplitudes of v ones. Assuming n = 3 — 5 and 
ri = Q, for these pulsars with riots = 50 from relation H28[) 
one get: 



e«J 1 = -0.76 - -0.68 (30) 

V nobs 

Therefore, the existence of such a remarkable group of pul- 
sars on the Uobs — Tch diagram argues in favour of signif- 
icantly high relative amplitudes of v variations - 50{t) in 
decomposition ((8]). 

Young pulsars in the area III are mostly the ones with 
negative e. Since the variational process is likely symmetric 
relative to evolutionary trend i/cv, young pulsars with posi- 
tive values of e should exist with Tch < Tch,ev and Uobs that 
is likely less than n. 

The pulsars of area I seem to be exactly such objects. 
They all have Uobs < 3 and their observed characteristic ages 
are most likely shifted towards lower values. Also, there is a 
deficit of pulsars in the area ifl (with Tch = 2 — 7 kyr) which 
may be a result of pulsars escaping from there to lower and 
higher Tch due to positive and negative e, correspondingly. 
So, the variations can be a reason for Uobs < 3 measured for 
the youngest pulsars known. We will return to the analysis of 
the youngest pulsars within our concept of cyclic variations 
in the Discussion. 

Finally, all remaining pulsars on the diagram - objects 
of the area IV - form two nearly symmetric branches with 
positive and negative Uobs- The variational amplitudes of i> 
of these pulsars are much larger than corresponding evolu- 
tionary values: I77I 3> 1. They are also affected by large am- 
plitudes of variations that may produce the oldest (with 
Tch > 10^ — 10* years) observed pulsars. Indeed, the largest 
pulsars' Tch seem to be physically unreasonable, as pulsars 



^ Moreover, these two pulsars with negative i> (PSR B1338-62 
and PSR B1610-50), while satisfying to the formal selection crite- 
ria described in Sec. 12.21 have relatively bad data sets - B1338-62 
shows significant glitch activity, and time spans for both of them 
are only about 200 days long. 



Note, however, that at least six "ordinary" pulsars without u 
measurement are known in this area. Their ages concentrate near 
Tch ~ 3 and 5 kyr as seen from the set of red vertical strokes in 
Figure m 
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Figure 7. The 71^),^ — t^^ diagram. It is similar to the i/ — u one shown in Figure [T] Correlation coefficients for positive and negative 
branches in logarithmic scale are 0.78 and 0.76 respectively. Red vertical strokes represent values of T^h of 1337 ordinary pulsars, displayed 
to show overall distribution of r^f^. The younger pulsars, particularly associated with SNRs, are located on the left part of the diagram 
(areas I, II and III). Their i) amplitudes are still quite small and are unable to produce negative observed nojj even in area III. However, 
the braking indices of pulsars there are obviously anomalous (noi,s < 50), and can not be explained by non-constant coefficient in spin- 
down law ^ as illustrated by blue dashed lines representing pulsar tracks with —K/K = 5 X 10* (a) and lO'^ (b) years, respectively. 
It suggests that these values are primarily affected by large u variations [A ~ 0.5 — 0.7). At the same time, the lack of objects in area 

II combined with small riabs of pulsars of area I means that observed braking indices and ages of the latter group are likely to be also 
affected by u variations. The observed Tch of the youngest pulsars seem to be less than evolutionary values, while Tch of pulsars of area 

III vice versa are greater. The amplitudes of u variations of pulsars of area IV are much greater than evolutionary values which produces 
two nearly symmetric branches of the diagram. 



shoul d cross their "death lines" in a few millions of years 
only (|Bhattacharva et"al] Il992l : iFaucher-Giguere Sz Kaspil 
l2006l ). 



4 PULSARS' SPIN-DOWN MODEL 

In the previous section we performed a statistical analysis of 
pulsars rotational parameters and concluded that observed 
pulsar spin-down is consistent with an idea of the existence 
of some cyclic variational process. Model-independent esti- 
mations show that this process operates on timescales as 
long as few hundreds or thousands of years, while rela- 
tive amplitudes of variations are comparable to secular, 
monotonous component of pulsar spin-down Vev 

Also, we found some evidence that v variations are sym- 
metric in respect to evolutionary term j>e„ . This fact may be 
useful to reveal the parameters of the evolutionary compo- 
nent, common for all pulsars. 

In the Section below we construct a two-component 
model of observed pulsars' rotational evolution which con- 
sists of evolutionary monotonous and cyclic components ac- 
cording to decomposition ([8]), and derive its parameters. 



4.1 Monotonous component of observed 
spin-down 

For the secular spin-down term of our model we will assume 
canonical power-law expression: 

Uev = ~Ku"^, (31) 

with K — const unique for each pulsar. On the other hand, 
we presume the constant evolutionary braking index n to 
be the same for all pulsars. The observed quantity noba and 
model parameter n are not equivalent and are connected 
with relation 1)28^ for each pulsar. 

Combining equations (O, ((8]) and H3ip . the evolutionary 
value of second frequency derivative for each pulsar with 
observed v and u may be written either as 

i^ev,i = nK^u^"~^ (32) 

or as 

Vev,2 = nK~ i^Y^T^) ^^^^ 

These expressions correspond to projections of an evo- 
lutionary trend onto either iz — u or iz — i' planes, respectively. 
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Figure 8. Distribution of relative deviations e simulated accord- 
ing to relation 1134 l l and distributions I I35II and II36I I witii (A) = 0.8 
and o-[A] = 0.2. Tlie probability P{e < —1) determines the prob- 
ability for pulsar to be observed with positive u. 



Moreover, when computed using actually measured timing 
parameters, these quantities may be treated as a statisti- 
cally independent ones, as they arise from combinations of 
different observables. 



4.2 Cyclic term of observed spin-down 

As a toy model, not necessarily exact but still represent- 
ing all the important properties, we assume below a simple 
harmonic form for an oscillating term of sum ([8} : 



e{t) = A cos 'p{t), 



(34) 



where A is a constant relative amplitude, and ip{t) is a phase 
of variations. 

Actual values of e for individual pulsars, as well as val- 
ues of coefficient K in the law (|3ip . are unknown. Therefore, 
we will analyse our subset in terms of its average, ensemble 
characteristics only. The parameters of e's distribution will 
be included in the model. 

Namely, to take into account physical diversity of indi- 
vidual pulsars, we assume A to be normally distributed with 
mean (A) and variance o-[A]: 



A^M{{A},a^[A]). 



(35) 



At the same time, the phase tp is obviously distributed uni- 
formly. 



<p ~ I7m(0,27r), 



(36) 



as it depends (as a modulus of ip = ipo + f{T, t) quantity 
over 27r) both on an unknown value of initial phase tpo and 
randomly selected moment of observations - pulsar' age t. 
It is also safe to assume the phases of individual pulsars 
to be uncorrelated both over the ensemble and with pulsar' 
parameters. 

Simulated e distribution for typical parameters of l^A) = 
0.8 and (j[A\ — 0.2 is shown in Figure [S] Values of e less 
than —1 correspond to the positive observed !> since !> = 
i'eu(l + e)- The probability to find a pulsar with v > 0: 
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Figure 9. Probability Ppos = Vie < —1) for a pulsar to be 
observed with positive !> for different combinations of parameters 
(A) and (j[A] of e distribution. 



for such a model depends on the parameters of the e distri- 
bution (A) and i7[A]. This dependence is shown in Figure |9l 

4.3 Evolution of an ensemble of pulsars 

Since we assume n to be the same for all pulsars, their trends 
i^evii'ev) on the i' — v diagram are characterized by individual 
coefficients K and initial frequencies vo only. These trends 
are straight lines in logarithmic scale (see Figure |4]) and do 
not intersect each other. Hence, the average behaviour of a 
pulsars' ensemble may be characterized by the line separat- 
ing the set of their evolutionary trajectories in half. Such 
trend corresponds to the median value of the coefficient 



K = M 



1 + e 



(38) 



where M[-] is a median over an ensemble and relations ([8} 
and (|26p are used. The value of if = K{n) can not be esti- 
mated from observables only since values of e for individual 
pulsars are unknown, but one can build a distribution of 
K assuming statistical properties of s within the ensemble. 
This distribution is a function of three parameters: 

k ^ K{{A),a[A],n) (39) 

It will be used below to build a maximum likelihood estima- 
tor for the parameters of the pulsars in the ensemble. 

4.4 Criterion for the estimation of model 
parameters 

If the v distribution is symmetric with respect to D^v, then 
the numbers of objects with observed less and greater than 
corresponding i>e„ should be nearly equal, both for the whole 
ensemble and for any subset defined by quantities uncorre- 
lated with phase (i.e. f or v). It provides a simple criterion 
for a model goodness-of-fit for an appropriately chosen i>eu,i 
(defined by equations (|32[) and (|33|l l the numbers of pul- 
sars observed with i> > Dev,i should simultaneously, for both 
i = 1, 2 be distributed binomially as 



Ppos = V{e < -1) 



(37) 



Bin{N,p=-), 



(40) 
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where A'' is the total number of pulsars in the ensemble. The 
same should obviously be valid for any sub-interval along f 
or i>: if the number of pulsars in fc-th sub-interval on i-th 
plane is A'^^fc, then the probability to have A'^^^ of Nik pulsars 
with iy > Uev is 



Pik — 



iN.k-N+y.-N+[ 



(41) 



Using this criterion and the dependence of A^^^ on the 
parameters of pulsars' spin-down model, an obvious maxi- 
mum likelihood estimator can be constructed to derive these 
parameters. 



constraints the parameters of e distribution, as probability 
Ppos for the pulsar to be observed with positive i> unlikely 
exceeds 1/N and corresponding probability 



p 



1 Ppos > 1 



(45) 



for the pulsar to be observed with negative ii therefore has 
to be large enough. 

The binomial probability P^^g to observe all A'' pulsars 
with negative values of v has therefore been added as a mul- 
tiplicative term to the likelihood function, and 

• finally we calculated logarithmic likelihood functions as 



4.5 Maximum likelihood estimator for the model 
parameters 

The estimation of TV^^ for any interval oiv oxv requires exact 
values of i>eu,j and, hence, known specific K and e for each 
pulsar. These values can not be derived directly from obser- 
vational data and therefore we have to hypothesize about 
them. 

We adopted the following routine for its estimation. 
Since no ad hoc information is available on e, we assumed 
it to be a random value distributed according to equations 
(|34p - (|36p . The specific K has been, in turn, replaced with 
median K value computed according to relation (|38|) with 
these e values in mind. This way, the corresponding esti- 
mates for the monotonous component of the second fre- 
quency derivatives for each pulsar from equations H32p - (|33p . 



= -AflogPn 



Y^^ogPu 



(46) 



r>2 2n-l 
U^v 1 = nK V 



and 



■ nK- - 



(42) 



(43) 



are in turn random quantities, uncorrelated to each other as 
they arise from combinations of different observables. They 
are, however, median estimates of its exact values due to 
functional form of equations, and therefore should roughly 
divide an ensemble of pulsars in half for adequately chosen 
model parameters - n, {A) and (t\A\. 

So, to construct maximum likelihood estimator we di- 
vided all pulsars into m = 6 intervals of A'^^ ~ 50 objects 
each (k — 1, 2, ...m) along i/ and v on the u ~ v and v — v 
planes, respectively. Then we calculated a likelihood func- 
tion as follows: 

• On the first step specific values of (A), (t\A\ and n have 
been chosen from the parameter space. 

• Then A'^ = 297 values of e have been simulated accord- 
ing to distributions (f55)) and (f5S)) . and 

• K(n, {j4),(j[j4]) have been determined. 

• The number of pulsars A'^^ with v > Vev.ih in fc*'' inter- 
val have been calculated and 

• used to compute binomial probabilities Pih defined (|4ip . 

• Then, the likelihood functions has been constructed: 



The indices i = 1 and 2 here mark (statistically independent) 
likelihoods for the planes v — v and v — v, respectively. 

We performed simulations of relative deviations e — 
Su/u for 297 pulsars 10^ times in a rott[f]for each set of given 
(A), cr[A] and n. In this way, the distributions (but not exact 
values) of d for every combination of parameters have been 
obtained. We have computed these distributions over a grid 
of values n, (A) and (j[A] in the intervals 0.5 — 6.5, — 1.2 
and — 0.4 with steps 0.05, 0.01 and 0.05, respectively. 



4.6 Likelihood for the null hypothesis 

Classical maximum likelihood technique presuppose a search 
of extremum of obtained £ within parameters space and 
then constraining the confidence regions assuming that dis- 
tribution of jC is known if the null hypothesis is valid. Typi- 
cally £ is a single- valued function of model parameters, and 
the p—% confidence region consists of a values of C for that 
P{Co < C) ^ p/100, where Co is a random quantity dis- 
tributed as a logarithmic likelihood within null hypothesis. 
The minimum of £ is obviously within this confidence re- 
gion. 

This method may be generalized in a straightforward 
way for the case when jC is not a deterministic function of 
parameters but a random quantity with known distribution. 

To do it, we simulated the distribution of Co by com- 
puting its values in a way similar to one used earlier for £.i, 
but using simulated according to binomial distribution 
with p = 1/2, and simulated uniformly distributed P„eg on 
the [(A'^ — 1)/N; 1] interval. Therefore, the null hypothesis 
corresponds to the exact symmetry of numbers of pulsars 
with i> greater and less than i>ev,i in each of m intervals 
described above. 

Then, using known distribution of £.,, we independently 
estimated the probabilities Pi = P{Co < Ci) and conflated 
them into t he resulting u n ited v alue V using t he method 
described bv lHill fc MM (|2010D andlHUil |2oi3): 



V1V2 



Plp2 + (1-Pl)(l-P2) 



(47) 



Wp. 



(44) 



This relation represents the cumulative probability that 
both likelihood functions have an extremum in the same 



We have also taken into account an absence of pulsars with 
positive v in an ensemble of TV = 297 pulsars. This fact 



° All simulations have been performed using the "Chebyshev" 
supercomputer of Moscow State University. 
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Figure 10. Plot of the confidence regions derived tiirough the maximum likelihood estimator. The number of slices for fixed o-[A] and n 
are shown. Blue, green and red contours represent 65, 95 and 99% bounds. The values of evolutionary braking index n, mean amplitude 
(A) and its dispersion cr[A], consistent with the symmetry of u variations, are at n ~ 2.5 — 4, {A) > 0.5 and cr[A] < 0.25. 



point of a model parameters space. Thus, the isocontours 
of V mark the bounds of a 100 ■ (1 — "P) percent confidence 
regions. 

4.7 Results of maximum likelihood estimation 

The maps of final conflated probabilities P are shown in 
the Fig. [To] as a collection of slices of the parameter space 
for flxed (j[A] (upper plots) or n (lower plots). The 99% 
confldence interval covers the range of n: 

2.5 < n < 4 (48) 

The corresponding range of (A) depends on the ac- 
cepted value of (7[A]. The decision about the appropriate 
amplitude and its dispersion can be made from analysis of 
Tiobs — Tch diagram provided above. Indeed, according to es- 
timation (|30p . typical relative deviations e are as high as 
0.7-^0.8, which corresponds to the a [A] ~ 0.1 in Figure [TU] 
So, we believe that the most adequate choice is: 

0.5 < (A) < 0.8, (49) 

and 

a[A] ~ 0.1 (50) 

Obtained values of n are in a good agreement with 
the n = 3 expected from simple magnetodipolar pul- 
sar braking theory with co nstant (or slowly evolving) K 
l|Manchester fc Tayloj[l977l ). 

The median evolutionary trend, assuming n — 3, rel- 
ative i> amplitude (A) = 0.65 and ct[A] — 0.1, is shown in 
Figure O 



4.8 Timescales of cyclic process 

The model of observed pulsars' spin-down does not include 
explicit relation between variational phase ip, pulsar age and 
timescale of variations T. To estimate the distribution of T, 
we will formulate it separately. Namely, let 

(f = (fio + Qt, (51) 

where Q — 2n/T = const. 

Since functions e{t) and ri{t) deflned earlier are not 
fully independent and related through equation (|13|) . the 
timescales (or frequencies) of a variational process for our 
model may be estimated directly: assuming secular spin- 
down law (|3ip . e — Acosip — Acos{Qt + ipo), e — 
~QAsin{Qt + ipo), for each pulsar: 

^_ u Uotsjl + Acosip) - n 
V A sin ' 

To get the distribution of all possible f2's for the ob- 
served set we simulated uniformly distributed phases (y3, as 
well as values of n and A ~ according to probability distri- 
butions within the confidence regions derived in Section [4. 71 
Negative values were rejected during this simulation as 
corresponding to physically impossible combinations of pa- 
rameters. Distribution of 57 simulated for a\A\ — 0.1 is 
shown in Figure [TT] Corresponding timescales are clustered 
inside a 5 — 500 thousands of years region. 

Also, the Figure [TT1 displays theoretical distribution for 
the timescales of NS precession caused by an anomalous 
braking torque (see Section [5.21 and equation (|57p ) for com- 
parison. These timescale values depend on a number of phys- 
ical parameters of NS, such as surface magnetic field, initial 
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Figure 11. The distribution of a timescales of long-term varia- 
tions, estimated according to the observed pulsars' frequency and 
its derivatives, for spin-down parameters values derived in Sec- 
tion l4.7l (see text for details). The timescales are clustered around 
5 — 500 thousands of years and are in a good agreement with sim- 
ilar distribution of NS precession timescales caused by anomalous 
braking torque (dotted line, see Section 15.211 . Dashed lines repre- 
sent timescales Tup and Ttyp according to estimations I I19I I and 
pijl . respectively, for !> = lO"^* s~^. 



period etc. We used distr ibutions of B and Po derived by 
iFaucher-Giguere fc KaspH (^006): log B and Po distributed 
normally with {logB[G]) = 12.65, aiogB = 0.55, (Po) = 0.3 
s and apg — 0.15 s. Masses and radii of NS have been taken 
distributed uniformly within 1.3 — 1.6 Mq and 8—12 km 
intervals, respectively. Initial magnetic dipole angles have 
been chosen randomly from [0; 7r/2] interval. 

Both obtained distributions are in a good agreement 
with each other and consistent with the idea of a long-term 
thousands-of-years variations of a pulsars' timing parame- 
ters. 



5.1 Ages of young pulsars revised 

Let's return to the diagram Uobs — Tch (Fig. 0- There are 
eight pulsars in area III with Uobs > 0, associated with young 
supernova remnants. Assuming their Uots are biased mostly 
due to !> variations, i.e. rj — 0, corresponding instant e can 
be derived as 

£ = V^/^obs ~ 1 (53) 
Therefore, their evolutionary characteristic ages 

Tch.ev = Tch(l + e) (54) 

and magnetic fields 

= -7^, (55) 
Vl + e 

where B is according to equation (|27|l . can be estimated, 
which are expected to be closer to the physical ones. 

We analysed data on three of these eight supernova rem- 
nants that have more or less independent and confident esti- 
mations of their ages, and summarized the results in Table[T] 
In general, corrected pulsars' ages become more consistent 
with that of corresponding SNRs in all of these cases. While 
the estimations of magnetic fields in all cases exceed 10^^ G. 
At the same time, if one assume initial periods Po for these 
pulsars, then the same corrections of characteristic ages can 
be provided by the proximity of the Po and the current pe- 
riod. 

A similar correction for the pulsars of area I is not so 
evident. Whether the measured Uobs < 3 of the youngest 
Crab-like pulsars are indeed biased from evolutionary n — 3 
due to the same variational process is not clear. However, 
according to equation (|29p . pulsars born with nonzero char- 
acteristic ages Po/2K ~ 10^ years. If so, then all pulsars in 
area I likely have Tch,ev greater than at least 2 thousands of 
years while their measured Tch < 2 kyr imply positive e. 

We believe that understanding of a possibility of sig- 
nificant biasing of standard estimators of pulsars' ages will 
induce new insights on the PSR-SNR connections, evolution 
and kinematics. 



5 DISCUSSION 

In this work we revise the problem of anomalous values of 
isolated radiopulsars braking indices. Since individual val- 
ues of i> are strongly dominated by non-monotonous com- 
ponent of the spin-down, we performed a statistical analysis 
of the ensemble of 297 objects. As a result, we concluded 
that pulsars' behaviour on the v — v and i/ — v diagrams 
suggest the existence of long-term cyclic mechanism affect- 
ing observed timing parameters. We constructed a semi- 
phenomenological model of observed pulsar spin-down and 
estimated it parameters. 

We found that the timescale of the v long-term varia- 
tions is likely close to few thousands of years - it follows both 
from simple model-independent estimation (Section l3.4|) and 
from accurate maximum likelihood analysis (Section 14. 8p . 

In general, our analysis inevitably suggests that ob- 
served V of radiopulsars vary with large enough relative am- 
plitude. As a result, basic pulsar parameters depending on 
V may be significantly over- or underestimated. We already 
said some words about this effect in Section [3.51 and below 
we'll continue its discussion. 



5.2 Long timescale precession of a neutron star? 

And finally we will discuss possible physical nature of a long- 
term variations of observed pulsars' spin-down suggested in 
this work. 

The unusual timescale of these variations - thousands of 
years - is much larger than the pulsar rotational period and, 
on the other hand, much shorter than its typical lifetime. At 
the same time, even the canonical magnetodipolar braking 
model seems to already include an essential mechanism for 
long-term variations. 

More than half a century ago iDeutscbl (Il955l ) derived 
the vacuum solution for the electromagnetic field of a per- 
fectly conducting, rigidly rotating spherical star. The full 
braking torque affecting such a star with magnetic moment 
pL and spin angular velocity may be described as a super- 
position of three orthogonal terms (in observer's frame): 

tl; = Q • (/2 X cj) + /? ■ ((/2 X oj) X a3) + 7 ■ (56) 

where \i = const. The second and third components of this 
braking torque are due to the magnetodipolar radiation in 
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Table 1. Corrections for characteristic ages and magnetic fields of several pulsars due to u variations. Here e is an estimated relative !> 
shift assuming evolutionary braking index n = 3 and 5u = 0, while Pq is an initial pulsar period that is necessary to provide the same 
correction. 
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far (wave) zone. B oth of them vary as (ujR /c)^, where R 
is the star radius l|Davis fc GoldsteinlHoTOI ). At the same 
time, first component of torque (|56p is due to radiation in 
near zone and varies as {ujR/cf . Since u)R/c « 10~* — 10~^ 
for a typical pulsar, the near-field torque is up to four orders 
of magnitude greater than far-field one. 

However, being very strong, this torque does not di- 
rectly affect pulsar spin-down rate and can not change mag- 
netic inclination angle x, as it is always perpendicular to 
both rotational and magnetic axes, in contrast to the far- 
field p art of the torque driving the evolution of pulsar' uj 
and X (|Davis fc GoldsteinllToTg ). 

This near-field torque induces an additional complex ro- 
tation of the star. Indeed, first term of torque H56p describes 
a precession of a neutron star around its magnetic axis with 
angular frequency = \a^\. Due to the strength of the near- 
field torque, the characteristic timescale of this precession is 
significantly shorter than the pulsars' typical lifetime. Pre- 
cisely, in a simple vacuum magnetodipol ar case, precession 
period remains constant and is equal to l|Good fc Ng|ll985l ) 

r = = (2.9x10' yr)-/45-/^-^o-B.:?2-^?6'-cos-'xo (57) 

Here Ui - pulsar's initial frequency normalized to 50 Hz; 
Bs,i2 - surface magnetic field normalized to 10^^ G; lis - 
NS moment of inertia normalized to 10*'' g • cm^ , and R - 
NS radius normalized to 10^ cm; xo ~ initial magnetic in- 
clination angle. The value oi Q. <^ uj depends strongly on 
the NS parameters, hence one can expect a wide range of 
precession periods, from hundreds to thousands years. The 
distribution of these periods for reasonable values of NS pa- 
rameters is shown in Figure [11] and is in a good agreement 
with estimations of variational process timescales derived in 
this work. 

Note, that similar long-term precession on a timescales 
of few thousand years can be also caused by the NS distor- 
tion along the m agnetic axis due to the strong magnetic field 
jGoldreichlllOTOI ). 

However, if we associate pulsar beam symmetry axis 
with NS magnetic moment then it is possible to show that 
such NS precession itself can not cause large observed de- 
viations of i) (and u) - moreover, it is unable to power a 
cyclic variations in observed NS spin-down rate at all. In- 
deed, the braking torque vector (|56p is dominated by the 
strong near-field term. Hence, assuming /? ~ and 7 ~ 0: 

to ^ a ■ {p. X u) — {Q, X Lu), (58) 

where is a constant precession frequency. At the same 



time, observed pulsar spin frequency is defined by rotation 
of p. However, this vector rotates around angular velocity 

fl — {uj X fl), (59) 

which also changes its direction relative to observer. There- 
fore, the value of observed spin frequency of jl is not equal 
to UJ and is not clear from equations (|58|l and (|59|l only. 

But, combining these two equations we found that u + 
aft — 0, which suggests introduction of the vector 

L = Lu + afl = 4- SI = const (60) 

which is useful for the description of fl rotation by a simpler 
equation similar to (I59|l : 

P^iLxp) (61) 

Since L — const, it is clear that a star precessing around 
its magnetic axis looks like a non- precessing star with a bit 
biased observed spin frequency 

L ~ (jj + i} cos X, (62) 

where X is a magnetic inclination angle. 

Therefore, even slow monotonous evolution of Q, u) and 
X are unable to introduce any irregularities in observed spin- 
down rate. The discussed type of precession is indiscernible 
for an observer. 

Therefore, if observed long-term spin-down variations 
are really caused by such precession, some geometrical or 
physical mechanism connecting rotation of around pL with 
modulatio n of x. etc should be sugg ested. 

Thus, iBarsukov fc TsyganI (I2OIOI ) considered an alter- 
ation of electric currents in the NS magnetosphere caused by 
such precession, whic h leads to modulation oiu). At the same 
time, iMelatosI l|2000l ) have shown that for non-spherical but 
biaxial or triaxial NS, near-field torque will strongly affect 
its rotation. If pulsar's rotational axis (and magnetic mo- 
ment) significantly deviates from the one of its main inertia 
axes, then a very large variations of x (and consequently 
of Co) occur. H owever, variations so large are not observed 
(|Melatosll2000l ). 

Generally, any change in the observed pulsar timing 
parameters are determined by the variations of the dip/dt, 
where ip is the dihedral angle between two planes - the plane 
of symmetry axis of a pulsar beam and the plane formed by 
rotational axis and direction to the observer. Detailed in- 
vestigations of observed complex timing behaviour should 
always take into account the geometry of pulsar emission. 
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Ultimately, whether the monotonous precession of NS 
rotation axis around the magnetic moment able to produce 
any non-monotonous peculiarities in observed v and /> is not 
clear, but at least this process operates on sufficiently long 
timescales of thousands of years. On the other hand, there 
are no obvious physical reason behind some hypothetical 
stochastic process with ultra-red spectrum, able to produce 
the same amplitude of pulsars' second derivatives variations. 
The regular nature of these variations seems phenomenolog- 
ically preferable. 



6 CONCLUSIONS 

In this work we analysed the properties of a set of 297 "ordi- 
nary" radio pulsars with a well detected, via a pulsar timing, 
second frequency derivative. As a result we demonstrate that 

• the long term spin-down of pulsars is well described 
by the superposition of a "true" monotonous and a long 
timescale non-monotonous cyclic term; 

• the subsets of these pulsars with positive and negative 
second derivatives are statistically different, both in shape 
and in number of objects; this effect is stronger for younger 
pulsars and vanishes for the older ones; 

• this difference reflects the existence of an evolutionary, 
secular positive trend ; the i> variations around this trend 
are nearly symmetric; 

• simple binomial arguments based on the assumption of 
symmetric variations around the evolutionary trend allow to 
construct reliable estimator for pulsars' long-term spin-down 
parameters; 

• pulsars secular spin down is consistent with a classical 
magnetodipolar braking with n ~ 3 

• the mean relative amplitude of the first frequency 
derivative {i>) variations is as large as (A) = 0.5 — 0.8 with 
variance a [A] ~ 0.1 

• the characteristic timescale of these variations is likely 
to be of several thousands years; 

• consequently, the observed values of pulsars' character- 
istic age are biased by the factor of 0.5 — 5 from the "true" 
secular value Tch,ev\ this fact naturally explains the discrep- 
ancy between real ages and characteristic ages of several 
pulsars associated with supernova remnants, the presence of 
pulsars with extremely large, up to ~ 10* years, character- 
istic ages, and low braking indices of the youngest pulsars. 

• there is at least one physical reason able to produce such 
variations in magnetodipolar model - a complex neutron 
star rotation relative its magnetic axis due to influence of 
the near-field part of magnetodipolar torque. 
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